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Large-amplitude behavior of the Grinfeld instability: a variational 
approach 
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Abstract. In previous work, we have performed amplitude expansions of the continuum equations for 
the Grinfeld instability and carried them to high orders. Nevertheless, the approach turned out to be 
restricted to relatively small amplitudes. In this article, we use a variational approach in terms of multi- 
cycloid curves instead. Besides its higher precision at given order, the method has the advantages of giving 
a transparent physical meaning to the appearance of cusp singularities and of not being restricted to 
interfaces representable as single-valued functions. Using a single cycloid as ansatz function, the entire 
calculation can be performed analytically, which gives a good qualitative overview of the system. Taking 
into account several but few cycloid modes, we obtain remarkably good quantitative agreement with 
previous numerical calculations. With a few more modes taken into consideration, we improve on the 
accuracy of those calculations. Our approach extends them to situations involving gravity effects. Results 
on the shape of steady-state solutions are presented at both large stresses and amplitudes. In addition, 
their stability is investigated. 

PACS. 47.20. Hw Morphological instability; phase changes - 05.70.Ln Nonequilibrium and irreversible 
thermodynamics - 46.25.-y Static elasticity - 81.10.Aj Theory and models of crystal growth; physics of 
crystal growth, crystal morphology and orientation 



• ^ 1 Introduction 

X 

When a nonhydrostatically strained solid has a surface, 
■ ■ ■ ' at which material can be redistributed by some appropri- 
ate transport mechanism, it may reduce its elastic energy 
via surface undulations. Intuitively, this should be clear: 
stresses are partially relieved in the maxima of corruga- 
tions and enhanced in their minima. The elastic energy 
density is therefore reduced in the maxima and increased 
in the minima, favoring growth of the former and deepen- 
ing of the latter. This mechanism is at the origin of a mor- 
phological instability leading to the formation of grooves 
with a relatively well-defined initial spacing under uniax- 
ial stress PI2 and, possibly, the evolution of islands, if 
the stress is biaxial |SE]- Pertinent transport processes 
are melting-crystallization for a solid in contact with its 
melt and surface diffusion for a sufficiently hot solid in 
vacuum. The latter case is relevant in epitaxial growth, 
where the lattice mismatch between different materials is 
the source of biaxial stress. 

The instability seems to first have been predicted by 
Asaro and Tiller 5 , but its universal nature was recog- 
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nized by Grinfeld 1^, hence it has often been referred to 
as Grinfeld or Asaro-Tiller-Grinfeld (ATG) instability. An 
unambiguous experimental demonstration of the instabil- 
ity was given by Torii and Balibar [S], using solid helium 
in contact with its superfluid. 

It should be emphasized that the surface undulation 
evolving as a consequence of the instability is not due to 
elastic deformation such as bending (as would be the case 
on application of a pressure to a long thin rod, leading 
to the Euler buckling instability). Instead, the instability 
materializes itself via mass transport and is independent 
of whether the stress is tensile or compressive. When the 
solid is in contact with its melt, the latter is a particle 
reservoir, rendering mass transport easy (and the dynam- 
ics is not governed by a conservation law) . When the solid 
is in contact with vacuum such as in the case of heteroepi- 
taxy, the instability takes place via surface diffusion in 
most cases, but may also be supported by other transport 
processes such as vacancy or impurity diffusion. In that 
case, mass conservation is important in the dynamics. For 
a pedagogical introduction into the subject of the ATG 
instability, we refer the reader to "Tl. 

There are a number of interesting questions concern- 
ing the instability. Since it produces crack-like patterns 
IS], does it constitute a generic route to fracture as hy- 
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pothesized in [HI or will plasticity in general lead to a 
restabilization? If one restricts oneself to linear elasticity, 
are there steady states beyond those found by Spencer and 
Meiron JOI? It has been shown that in directional solidi- 
fication stable steady state patterns are realizable [TItTT] . 
For the pure Grinfeld instability, this appears to be impos- 
sible in extended systems. Further questions concern the 
nature of dynamical states, when there is no steady state. 
Coarsening has been found to be a generic behaviour |12[ 
Il[{|ll4j , but more detailed investigations on large-scale sys- 
tems would be desirable, to determine the precise form of 
the pertinent power laws. 

Numerical simulations of a solid undergoing the Grin- 
feld instability (S^iQ-lO, have the awkward tendency of pro- 
ducing cusp singularities in finite time. The investigation 
of Spencer and Meiron lOT has shown that these singu- 
larities are not an artifact of the numerics but intrinsic to 
the continuum model describing the system, under the as- 
sumption of linear elasticity (and in the limit of negligible 
sound propagation effects). What they found was a steady 
branch of solutions in a certain range of wavelengths, cor- 
responding to very small sinusoidal shapes near the onset 
of the branch and approaching a cycloid-like cuspy shape 
near its termination. 

Such a result might have been anticipated on the ba- 
sis of the analytic work by Chiu and Gao who per- 
formed a detailed calculation of the stress state under a 
cycloid-shaped surface using the Goursat function scheme 
proposed by Mufichelischwili ^B). All of these numerical 
studies considered two-dimensional systems where the in- 
terface is described by a curve. Whereas we have treated 
three-dimensional systems in 0] within a weakly nonlin- 
ear approach, we will restrict ourselves to two dimensions 
here but go well beyond the regime of validity of weakly 
nonlinear theory. 

Chiu and Gao find that for a certain range of wavenum- 
bers a fully cusped cycloid constitutes an energetically 
more stable configuration than a flat surface. In section 
12 we will show that a variational calculation using cy- 
cloids as ansatz functions gives a rather good approxima- 
tion of steady state solutions of arbitrary amplitude, some 
of which were discussed in jlOj . 

Moreover, we are able to draw conclusions on the large 
amplitude behavior even for strong gravity or a large den- 
sity difference between the solid and nonsolid phases (liq- 
uid or vacuum), as we show in section |21 Evidence for 
these states has already been found in |17| . 

In section 0] we present a generalization of this idea. 
Employing a special system of (not necessarily univalent) 
functions called multi-cycloids we analytically recover the 
numerical results for the mean square amplitude to ex- 
cellent accuracy already at third order. At higher order, 
we get more precise results with less numerical effort than 
Ref. [ini. 

Finally, we give some conclusions as to the physical in- 
terpretation of our results and suggest how to verify them 
experimentally or by a full numerical computation. 



2 The mono-cycloid approximation 
2.1 Cycloids 

We wish to use cycloids and more general curves deriving 
from cycloids to model the steady-state surface pattern of 
a two-dimensional solid after it has undergone the Grinfeld 
instability. This is of course motivated by the fact that 
cycloids have been shown by Chiu and Gao ^S] to be very 
efficient in reducing the elastic energy, and hence the final 
steady state, if any, should be close to a cycloid shape. 

Cycloids belong to the more generic class of trochoids, 
curves defined as the trace of a point fixed on a circle 
rolling along some prescribed line. A cycloid is the curve 
traced out by a point on the circumference of a circle as the 
latter rolls along a straight line. When we put the point 
inside or outside the circle instead, we obtain a curtate or 
a prolate cycloid, respectively. The parametric represen- 
tation of a cycloid can be given in a compact manner by 
a complex generating function 

C(C)=?-*fe-*'«, (1) 

Herein, i is the imaginary unit, /c is a wavenumber and p 
is a dimensionless amplitude-like parameter. The cycloid 
is obtained by taking the real and imaginary part as the 
X and y coordinates, respectively. 





9(C) 
-1 



-2 




Fig. 1. Cycloids, as in equation ^ with k — 1, plotted in the 
range ^ — {—it, tt); p = . . . 1 from top to bottom. Curves have 
been shifted in order to avoid overlapping. 

Taking p = 1 leads to the representation of the clas- 
sical, i.e., cusped, cycloid. Choice of a plus instead of the 
minus sign in equation ^ would just shift the minimum 
from ^ = 0to^ = 7r/fc, while a plus sign in the exponent 
would lead to a surface with the cusps pointing upward. 
The latter case is of no relevance for the Grinfeld insta- 
bility, because a fully relaxed upward cusp would immedi- 
ately shrink under any perturbation in order to decrease 
the surface energy. 

For p > 1, the cycloid becomes self- intersecting and 
does not represent a physical state anymore. Note, how- 
ever, that if we "superimpose" several cycloids^ as we will 

^ This is not a superposition in the standard sense. We add 
up partial representations of the x and y coordinates, so the 
whole curve is not a simple sum, see Eq. I|29|l . 
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do in section 0] the possibility of the x coordinate to 
vary nonmonotonously allows the representation of pat- 
terns with overhangs that do not self-intersect. We will 
not discuss this feature in detail here but report on those 
patterns in a different article. 

The parametric representation of the cycloid is 



xiO =3?[C(0]=e-|sin(fcO, 



■ cos (fc^) 



(2a) 
(2b) 



Assuming the surface of a solid having undergone the 
Grinfeld instability to be described by this shape, we shift 
the cycloid position by its mean value 



k 

2^ 



El 

2k ' 



(3) 



(To = '^xx — <^zz is the first normal stress difference or, in 
more physical terms, the excess stress applied in the x di- 
rection, to produce a uniaxially strained solid, whereas 
E and v are elastic constants describing an elastically 
isotropic material, viz. Young's modulus and the Poisson 
number. 

Physically, the Griffith length describes the competi- 
tion between surface energy and elastic energy. It is used 
predominantly in the theory of crack propagation. Cracks 
larger than this length relieve more elastic energy when 
growing than they produce surface energy, while cracks 
shorter than it can reduce energy only by shrinking. There- 
fore, this length scale represents a nucleation size for crack 
generation. 

When gravity is considered, another length becomes 
important, which is 



in order to keep the average interface position at ?/ = 0. 
Herein, the dot denotes differentiation with respect to ^. 
Hence, from now on we use 



y(0 



(4) 



Moreover, if we want to compare results with both our 
amplitude calculation J7] and the work of Spencer and 
Meiron ^U], we additionally have to know the mean square 
amplitude of the cycloid. As in ^7], we denote this mean 
square amplitude by a ; its calculation yields 



k 

2^ 



- 1/2 


Pv/2 - 




2k 



(5) 



where the absolute value ensures a nonnegative mean square 
amplitude, no matter what the sign of p and of k. 



2.2 Scaling 

Our basic model is an isotropic solid obeying the laws 
of linear elasticity (i.e., the Lame equations) with a sur- 
face on which shear stresses vanish while the normal stress 
component is equal to the negative pressure in the liquid 
or zero. That is, we neglect capillary overpressure due to a 
curved interface, which is known to be a small cross effect 
U]. Moreover, we neglect the body force effect of gravity 
in the elastic equations, also known to be small. 

The energy of the solid then consists of three contribu- 
tions, its elastic energy, total surface energy, and potential 
energy in the gravity field. 

It is known that if the latter force, gravity, is com- 
pletely neglected, the equations of motion of the Grinfeld 
instability can be made parameter free. This is achieved 
by referring all lengths to a length scale li, essentially the 
Griffith length, given by: 



^1 = 



7 
2wo 



(6) 



where 7 is the surface tension and wq = crQ(l — i''^)/2E 
the elastic energy density of the prestressed planar state. 



Wq 



(7) 



Herein, Ap is the density difference between the solid and 
the second phase, g the gravitational acceleration. This 
gravity length describes competition between elastic en- 
ergy and potential energy in a gravitational field. Due to 
the density difference between the two phases, the system 
can gain potential energy, if the phase with the larger den- 
sity, usually the solid, shifts its center of mass downward. 
As a consequence, if a solid immersed in and in equilib- 
rium with its melt is submitted to a uniaxial stress, it will 
first start to melt, because it is now out of equilibrium; 
but then its center of mass shifts downward, hence the 
potential energy is decreased and a new equilibrium state 
may be reached. This happens whenever the applied stress 
difference is below the instability threshold. The solid sur- 
face melts back by a certain height, and this height change 
is exactly given by 12- 

The only parameter of the nondimensionalized equa- 
tions is the ratio of these length scales: 



/ - '1 
'12 •— Trr 



(8) 



In all considerations that follow, we have carried out a 
formal transformation x ~^ x/h, y y/h, k — > kli, ^ —^ 
^/h, and energies and their variations have been divided 
by a common prefactor 7. 



2.3 The cycloid approximation for the no-gravity case 

In treating the cycloid approximation, we choose an ap- 
proach that is essentially variational in nature. Therefore, 
we need not compute the energy itself but only its vari- 
ation. The variation of strain energy due to a configura- 
tional variation Sx may be written as 



J w{s)n6x ds. 



(9) 



where w(s) is the energy density at the surface, n is the 
normal vector and s denotes the arclength along the sur- 
face (as we are dealing with a two-dimensional system). 
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Using Eqs. H2a(l and we can calculate an approxima- 
tion to 6Ec, allowing only the parameter p to vary (instead 
of taking the full variational derivative which would take 
the result out of the space of cycloidal shapes) 



ox = — — e 



dy{0 



Sp 



\ dp "'^ ' dp 
= - k^^ [sin(fc^)ej; + (cos(fc^) + p) ey] Sp , 

and we have the relation 

nds = {-y (^) e^; -I- i {£,) By) d$, 

= [(-psin(fc^)) + (1 - pcos(fc$)) By] d£_ 



(10) 



(11) 



The calculation of the strain energy density can be per- 
formed for the general case of a multi-cycloid. Essentially 
the same calculation has been carried out before by Yu 
and Suo ^Hl , who used it to model groove-to-crack evolu- 
tion in ceramics, a context quite different from ours. Since 
the notations used by the two groups are pretty different 
and a direct translation would be tedious, we present the 
important steps of our calculation (done independently 
and based on rather than ITF) in appendix At 
this point we only need the energy density for the mono- 
cycloidal surface: 



1 



w(s) = — 

2(i + p2_2pcos(fce)) 



2 • 



(12) 



When comparing with the results of 15 , one finds the only 
difference (up to prefactors) in the different sign of the 
cosine function, which results from their different ansatz 
of the generating function, corresponding to a simple shift 
of the argument of the cosine. The integration yields a 
surprisingly simple result: 



dp 



^ ^«;(s)cos(fcC)(p2-l)de--^, (13) 



where we have taken fc > (otherwise the expression on 
the right-hand side would have to be multiplied by the 
sign of fc). 

It remains to be noticed that the integration can be 
done analytically in the mono-cycloid case, yielding 



irp 
'"fc2" 



(14) 



The surface energy in our scalings simply is the dif- 
ference of the arc lengths, and its derivative is calculated 
straightforwardly : 



E. 



i(0 +2/(0 d^-l 



k \l + pj 



dEs 
dp 



P-I^f2^p 



(15a) 



(15b) 



(16) 



Herein, K and E are complete elliptic integrals of the first 
and second kind, respectively, defined by 



K{u) 



E{u) 



V2 



dx 



2 ,• 2 , 



yl — u^ sm" X 

V2 



\u\ < 1) , (17a) 



dxVl-u^ ain^x, (|u| < 1) . (17b) 



The result l(TB|l simplifies to d^^/dp = ^ju in the fully 
cusped limit, i.e., for p = \. For later simplifications we 
set 



^{p) ■■= - ({p-I)k(^) +{p+l)E 

TT V yi + p/ yi + p 



2Vp 



leading to 



dEs ^ N{p)Ti 
dp kp 



(18) 




Fig. 2. The function N{p) and its derivative 



In Fig. 121 we plot the function N{p) and its derivative 
dN{p)/dp = 4pK(2^/(l-|-p))/[7r(l + p)]. First, we note 
that N{p) is monotonously increasing from A^(0) = to 
A^(l) = Ytt. Second, the derivative diverges logarithmi- 
cally near p = 1. At p — 0, N{p) is regular, and the first 
few terms of its Taylor series are given by 



Nip) 



8^ 



This expansion will become useful later in the discussion 
of the type of bifurcation at the instability threshold. 

To obtain the energy minimum, we simply have to 
solve 



dE, 
dp 



dE, 
dp 



2tt 

T 



P Njp) 
k 2p 



= 0, 



(20) 



The complete branch is obtained by solving Eq. 1)20(1 for 
k instead of p, taking into account that p is running from 
zero to one: 



Nip) • 



(21) 



The solution, additionally converted to a via equation (01 
for easier comparison, is shown in figure |31 At the termi- 
nation point of the numerical solution |10| , the a value of 
the mono-cycloid approximation is about 10% smaller. 
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T 1 r 

1.6 1.7 1.8 1.9 2.0 

Fig. 3. The solution of the mono-cycloid model (irl^ in compar- 
ison with the digitized and scaled solution branch of Spencer 
and Meiron (s). 



From (|21|l . we conclude that the analytical termination 
point, given by the cusp limit p=l, is located at k = ^/2, 
which is somewhat different from the termination point of 
Spencer and Meiron 10 at about 1.74 (their scalings im- 
ply double wave number and half amplitude with respect 
to ours, so they gave the termination point at fc = 3.48). 
We shall see later that the wavenumber of the true cusp is 
closer to the monocycloid result than to the numeric value. 
This is due to the fact that a very small change of tip ra- 
dius in the minimum of the pattern leads to a relatively 
large change of the wave number. So the Spencer-Meiron 
result is accurate for the part of the solution branch that 
it reproduces, but it misses a considerable piece of the 
branch. 

From Eq. Q we find using fctcrm = 



(22) 



An important question is the stability of our solutions. 
Again, only a simple calculation needs to be performed: 
Take k from equation (|21|l . and insert it into the second 
derivative of the total energy E — + . This leads to 



k d^E 



N(p) , N'{p) Nip) , N'ip) 



-E 



2p2 



2p 



1 



il + p)p' 



p2 2p 

Km: 



(23) 



which is positive for any p e (0, 1). Hence the complete 
solution branch is stable up to the singularity, a result that 
agrees with that of Spencer and Meiron JlO^. Of course, 
with our method statements about stability can be made 
only concerning the restricted set of functions used in the 
variational ansatz (which depends on just one parameter 
here) . 

Beyond p — 1, Eq. flU^ still gives us a stationary point 
of the energy (even a minimum for p not too far above 
1), but the corresponding cycloid self- intersects, hence the 
solution is unphysical. Therefore, the variational ansatz 
provides a transparent analytic explanation of the fact 
that the solution branch terminates. 



3 Including gravity 

The next natural step is the incorporation of gravity into 
the model. Again, the calculation of the corresponding 
energy contribution is fairly easy, because up to a prefactor 
it is nothing but the square of the mean square amplitude. 

E.^'flmy((f^='-^!^^ (24) 

Consequently, the derivative is 



dE„ ttIuP (1 - p^ 



dp k^ ' 

and for the generalized model we replace Eq. (|20|l by 
d 



dp 



(E, + E, + Eg) = 0. 



(25) 



(26) 



Again it turns out that the solution can be written down 
exactly if we express k through p. Now we have two solu- 
tions which are both meaningful: 



P 'l±Jl~N{p)- 



N{p) 



p- 



tl2 



(27) 



As in the no-gravity case, we construct solution bran- 
ches by fixing I12 and calculating {k,p) pairs which may 
then be converted via Eq. (j^l into (fc,a) pairs. 

Only the solutions for I12 G [0, 1] cover the whole range 
p = . . . 1, while in the region I12 > 1, i.e., below the crit- 
ical point, the system is linearly stable and hence we find 
no solutions close to zero. In these cases it is necessary to 
first calculate the minimum possible value of p by requir- 
ing the radicand in Eq. (|27|) to be zero. 

Clearly, the fact that finite-amplitude solutions exist 
at subcritical values of I12 is already indicative of the sub- 
criticality of the bifurcation at the threshold, a result first 
obtained by Nozieres ^Hl- Let us discuss the vicinity of 
the critical point in some more detail. The neutral mode 
emerging at that point has of course p = 0. So we should 
expand the energy or its derivative for small p to obtain 
the solution behavior at the bifurcation. Using Eq. Ijl9|) . 
we find for the derivative of the total energy: 



k dE 
2^~dp 



1 

"k 



1 

16 



in 

2fc2 



2fc^'^ 
< 3 



128 



P 



25 
2048' 



(28) 



At the bifurcation point, the linear term vanishes, because 
k = 1 and I12 — 1. The third-order term is — ^/le p^, i.e., it 
is negative, whereas higher-order terms are positive. This 
is to be contrasted with the Nozieres calculation and its ex- 
tension by ourselves in JTj , showing that all the calculable 
coefficients of an amplitude expansion in terms of Fourier 
modes are negative. From this phenomenon we concluded 
that there is no restabilization at finite amplitude. At first 
sight, the situation seems to be different here, due to the 
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positivity of higher-order coefBcients. But in fact, it is not, 
because the maximum meaningful amplitude is p = 1, and 
it is easy to see that for this value the third-order coeffi- 
cient remains dominant. As p ^ 1, the energy derivative 
tends to the negative value 4 — 2tt. In a sense, this consid- 
eration shows a little more than the calculation in terms 
of Fourier modes: restabilization of the structure would be 
possible at p > 1, but this is an unphysical situation. 

Figure 01 shows an example for the profile of a steady 
state-solution taken from the critical branch {I12 = 1). The 
solution shows the typical behavior predicted by Nozieres, 
i.e., flat cell tips and sharp grooves. It is however unstable, 
as shall be seen from the discussion below. 



At even higher values of I12, we are in the subcritical 
range where it takes some energy to modulate the surface 
in order to overcome the energy barrier to let the insta- 
bility emerge. In Fig. we have shown solutions up to 
li2 = 2. The solution branches converge towards a = i/2fe 
as li2 is increased, which is the transformation of the cusp 
limit p = 1 (see equation ©). The common endpoint of 
all curves {'^/2, i/tt) is marked with a circle. 

Again, as in the no-gravity case we have to check the 
stability of the solutions by inserting the corresponding 
solution (|27|) into the second derivative of the energy. First 
we note that all branches left of A; = ^^72 in Fig. [3 start off 
unstably. But there is a range of stability which we discuss 
with the help of Fig. rather than|S| 



0.2 
0.0 

-0.2 
-0.4 




0.00 2.42 4.83 7.25 ^, 9.67 



Fig. 4. Two periods of a sample solution at A: = l^S^x? = 0.25 
(which corresponds to p « 0.49). 



A number of solution branches for different values of 
1x2 is shown in Fig.|31 In this plot, we first of all see how 
the range of unstable wave numbers, starting with the 
interval (0,2) becomes narrower with increasing Z12, i.e., 
higher gravity or lower prestress. At I12 = 1, only a single 
k value has a non-negative growth rate; this is the critical 
case. 




Fig. 5. Solution branches of equation 12611 for differ&it values 
of Z12. The rightmost curve, beginning at (fc, a) = (2.0,0.0), 
is the no-gravity solution shown in Fig. |3| Arrows and small 
numbers denote Z12 values corresponding to the curves. The 
upper thick line is the limiting curve for large Z12. The common 
termination point of all branches is marked by a circle. For 
more details of this region see Fig. El 




Fig. 6. Solution branches of equation 12611 in a certain k- 
range. The region bounded by the no-gravity solution and the 
curve (thick line) from the bullet symbol up to the cusp point 
contains the stable solutions. 

In this figure, we have included the range of stable 
amplitude solutions. The largest I12 value where small so- 
lutions still emerge stably (•) is I12 = 32/8i, correspond- 
ing to fc = This should be compared with the re- 
sult of jl7|. where the tricritical point is shown to be 
located at fc = 13/3 — ^^/a ~ 1.82, corresponding to 

= 20V57/9 _ 148/9 « 0.33. Because the result from 
the amplitude equations is exact at infinitesimal ampli- 
tudes, the variational ansatz overestimates the range of 
supercritical bifurcation at lowest order. Of course, such 
a comparison is a bit unfair towards the variational ap- 
proach, as we have only a single parameter available there, 
whereas the lowest-order nonlinear mode expansion con- 
tains already two amplitudes. As soon as we take the 
multi-cycloid ansatz, discussed below, to second order, we 
obtain the exact position of the tricritical point within 
this ansatz as well. This is simply due to the fact that the 
n-th summand in the multi-cycloid expression (Eq. I29|l 
does not contain Fourier modes lower than n, hence all 
contributions of order 2 must be present for N — 2 (but 
some are already present for iV = 1). 

An interesting fact is the bending of all solution bran- 
ches into a stable region before running into the cusp, 
a feature not obtained within the amplitude equation ap- 
proach. Actually, the stability changes at the points where 
the slope turns negative (passing from +00 to —00). This 
can be interpreted as a hint for the existence of stable so- 
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lution branches at high gravity. These stable solutions do 
not, at least in the mono-cycloid approximation, emerge 
at the upper marginal wave number 1 + Tl2 predicted 
by the linear stability analysis but instead in the (upper) 
vicinity of fc = ^2 at nonzero amplitude. Whether this is 
really true, will be checked in the following chapters by 
the implementation of the multi-cycloid ansatz. 



4 The multi-cycloid approximation 

In seeking a generalization of the cycloid ansatz, we were 
guided by two principles. First, the generalization should 
reduce to the cycloid, when all but one of the parameters 
became small. Second, it should correspond to a boundary 
curve that allows us the analytic solution of the elastic 
problem by conformal mapping. 

Given these conditions, the multi-cycloid model is a 
natural generalization of the cycloid one (and of a two- 
cycloid ansatz already considered in |14|'): 



c(e) = e- 



N 

^ nk 

n=l 



(29) 



Herein, N is the number of "cycloid modes" taken into 
account. The denominator n in this equation which could 
also have been included into the definition of the pn has 
been explicitly written in order to be able to express the 
generalized cusp condition in a compact way. Real and 
imaginary parts read: 



N 



y(C) 



N 

^ ^ ^^^^^ 



(30a) 
(30b) 



Again we shift the interface by its mean value to set the 
average interface position equal to zero: 



m 



2n I 2k ^-^ n 

n=l 



(31) 



and we have to correct Eq. (|30bp as follows: 



N 



(32) 



Moreover, we will need the mean square amplitude again. 
The result is 



N 



\ n=l 



2n2 



N ^2 

2 ^ j 



JV-r 



PnPjPn+j 

fr{ Jin+j) 



(33) 

We assume the sequence of the p„ to decrease sufficiently 
fast so the sum of the absolute values of the p„ does not 



exceed one, a condition that is sufficient to avoid self- 
crossings of the curve given by H29|) . Then cusps, if they 
exist, can appear only at 



Ccusp k = 27rn, neN . 



(34) 



(They are characterized by C'(0 = 0-) The radius of cur- 
vature is given as [201 



2\2 



m 2/(0 



(35) 



and it takes its minimum value at C = i^cusp- Therefore, 



N 

-1+ E pn 



N 

k E 



(36) 



and the cusp condition reads 



N 



(37) 



Now the derivatives of the relevant energies have to 
be calculated. Note that the integration for the energies 
themselves may not be carried out analytically for general 
multi-cycloids. But the question of stability of the solution 
can be answered, because first and second derivative can 
be given explicitely. 

The crucial point is the calculation of the elastic energy 
density w{s) at the surface, which is done in appendix IXI 
Supposing w{s) is given (see Eq. (|59|l ^. we have to modify 
Eqs. (|10ll and as follows (we write C„ and Sn instead 
of cos(nfc^) and sin(nfc^)): 



N 



5pn 



Opn Opn J 

1 ^ 

6p 



nk 

n=l 

nds = (-y (0 e^; -I- i {£,) By) 



Consequently, we get according to Eq. © 

Hs) [Sr^m~{Cn+Pn)±md^ 



dp,. 



1 

nk 



(38) 
(39) 

(40) 



The other terms are simpler again. Gravitational en- 
ergy is /12/2 times the square of the mean square ampli- 
tude lO 

E,^^fa\ (41) 

if we divide the integral by the wavelength, as we will do 
in all energy expressions from now on, i.e., rather than 



8 



Peter Kohlert et al.: Large-amplitude behavior of the Grinfeld instabihty: a variational approach 



integrals over a periodicity unit we consider averages. We 
obtain 



dPn 



ll2_ 
2fc2 




1 J't^'i 



J\ 



(42) 



and the surface tension is again represented by the dif- 
ference of the arc lengths (compare Eq. UlSal) ). hence its 
derivative reads after some simplification 



dpn 



Pn Cn+ J2 ClPn-l + ClPn+l 



1=1 



1=1 



N 



N-1 



N-l 



fiP, = J2p' - 2^') + 2 E E P^P^+^ (43) 



1=1 



1=1 



This integral can only be solved analytically for the case 
= 1; the analytical result has been used in the mono- 
cycloid model. 

Equipped with the terms (|^ . (|^ and we can 
now numerically solve the system 



OPn 



n = 1 . 



.N . 



(44) 



for the set {pi . . . pn}, given a certain prescribed k value. 
Some technical details of the solution method are described 
in appendix IbI 

We have carried out the calculations for a range of N 
up to 50. Figure 13 shows how fast the effective amplitude 
branches converge to the numerical result from |1(J| . 

A means to assess the closeness of a solution branch to 
its cusp termination is to determine the radius of curva- 
ture in the minimum of the grooves which will approach 
zero near the cusp. Here we note that the iV = 3 approx- 
imation is insufficient for a description of the tip radius 
(Fig.|SJ), in contrast with the effective amplitude (Fig.O, 
which is already well approximated by three modes for 
most of the branch. (We have omitted N = 2 here be- 
cause it is much worse.) 

As can be seen, the solutions agree with the Spencer 
solution in the range they cover. Yet, the extension to the 
range of lower wavenumbers is rather sensitive to the num- 
ber of included modes. A prolongation of the curves in Fig. 
IHl suggests a termination point slightly left of fc = 1.55, 
which is less than as proposed by the monocycloid 
model. As with increasing N the termination point moves 
to slightly higher values of fc, it is tempting to speculate 
that the exact termination point is at fc = 7r/2 indeed. In 
any case, our method allows to reach radii of curvature 
that are three orders of magnitude smaller than the min- 
imum value found by Spencer and Meiron, and it does so 
apparently with less numerical effort. 

The incorporation of gravity has been carried out up 
to li2 = 1. Figure El gives the solution branches, unsta- 
ble curves left of fc = 1 are not shown. We note that the 




Fig. 7. Comparison of the A*' = 2, 3 and 10 multi-cycloid 
models in the no-gravity case. The dotted line is the N = 2 
approximation, which is insufficient for our purpose: it does 
not reach the cusp at all but turns into a bag-like morphology 
instead (i.e., it develops overhangs). The dashed line is A'' = 3, 
being already in good agreement with jlUj . and the solid line 
is the A'^ = 10 example. All solutions with A'' > 4 look the 
same. Curves are terminating slightly before the cusp emerges, 
which will be further clarified in fig|Sl The crosses represent 
the numerical solution from |10| . 



10' 



10" 



10 ■ 



10' 



10" 




1.5 



Fig. 8. Radii of curvature for N up to 50, explanations see 
text. Again, the crosses show the digitized and scaled radii of 
curvature of the Spencer solution. 



common cusp point, an important feature of Fig. [HI disap- 
pears. As in Fig. [3 the curves terminate before the cusp 
is actually reached. 

Stability of the solutions is checked via computation 
of the determinant of the matrix d^E/dpndpm and its 
principal minors. If all of these are positive, the matrix is 
positive definite, the energy has a minimum and the solu- 
tion is stable. As it turns out, the determinant itself gets 
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positive only after its minors, so it would in our case be 
sufficient to check the sign of the determinant alone. The 
range of stable solutions is displayed as the grayed area 
in Fig. El We conclude that the stability behavior is qual- 
itatively correctly described by the mono-cycloid model 
already. However, stability sets in at smaller amplitudes 
for k values below 1.7. Hence, the multi-cycloid modes act 
stabilizing at larger amplitudes and destabilizing near the 
tricritical point, which is shifted to larger k by inclusion 
of the second mode. 



0.5 




1 1.2 1.4 1.6 1^ 1.8 2 



Fig. 9. Solution branches based on a = 30 multi-cycloid 
approximation. Numbers denote the value of Z12. The shaded 
area indicates the stable part of the solution manifold. Small 
sub-figures represent example morphologies. 

The subfigures of |^ show that up to the different dis- 
tances to the cusp, which make the curves more or less 
" sharp" , there is no indication from the shape itself whether 
it is stable or not. 



5 Summary 

To conclude, we have presented a variational approach to 
the calculation of steady states for the Grinfeld instabil- 
ity. Taking into account a single mode we already obtain 
a very nice qualitative description of the system behavior 
including the approach to a cusped state. The wavenum- 
ber for the cusp appearance is already more accurate with 
a single mode than in the article of Spencer and Mei- 
ron, while the amplitude is pretty far off the true result 
(by about the same amount as the amplitude obtained by 
Spencer and Meiron). 

Nevertheless, this single-mode approximation has the 
virtue of great transparency. That a cusp singularity ap- 
pears is rendered understandable: the system simply draws 
near a state where further minimization of the elastic en- 
ergy would require the interface to self-intersect. 

A few words may be in order concerning the limits of 
validity of our approach. The nature of our calculation 



is variational, which means that it will overestimate the 
energy of the system. Moreover, the minima of the vari- 
ational energy will not lie exactly at the same positions 
in parameter space as those of the true energy. As we in- 
crease the number of modes, we will get closer to the true 
result, and if our function system were complete, we could 
be certain of full convergence of the variational results to 
the correct answer. We have no formal proof of the com- 
pleteness of the system of multi-cycloids but note that as 
a function of ^, the systems used for the representation of 
the abscissa and the ordinate of the curve are complete in 
the spaces of odd and even functions, respectively. Com- 
pleteness is difficult to prove because of the correlation 
between the coefficients describing the abscissa and the 
ordinate. However, we suspect that for all practical pur- 
poses of representing curves that resemble a cycloid as 
closely as do the numerically obtained solutions, our func- 
tion system can be considered complete. 

In the full numerical computation |10| . with which we 
compared our results, the discretization of abscissae is 
given by a formula akin to H2a|l with an equidistant dis- 
tribution of the parameter ^ and the interface position 
is given as a superposition of cosine modes in the same 
parameter. Hence, numerical convergence relies on the as- 
sumption of completeness of a function system derived 
from Fourier modes by a stretching in the x coordinate. 

As long as the higher modes have small enough ampli- 
tudes, the two approaches should give equivalent results. 
Since we solve the elastic problem essentially analytically 
and for a continuous interface, not a discrete one, we reach 
the same accuracy as the numerics with fewer modes. 

Note that the cusp singularity is not an inherent re- 
striction to the method, as the function system is chosen 
such that it can represent one (or several) cusps. How- 
ever, when a cusp appears, quantities such as the elastic 
energy density diverge there. This means that the numer- 
ical solution of the nonlinear system of equations (|44(l for 
the variational parameters will run into problems, hence 
the cusp cannot be reached exactly in this final numerical 
step. 

Even the one-mode approximation suggests the exis- 
tence of stable large-amplitude steady states in the pres- 
ence of gravity, as is demonstrated by Fig. |H| Taking into 
account more modes, we obtain a quantitatively satisfac- 
tory description of the numerical Spencer-Meiron branch, 
conveying some confidence that the new branches with 
gravity are equally well described by this approach. The 
stability domain suggested by the one-mode picture is 
roughly confirmed in the three-mode representation (see 
Fig.inj. As gravity is increased, there are no small- amplitu- 
de stable solutions anymore. At first sight, this might seem 
counterintuitive: why should gravity, a stabilizing effect, 
destroy the stability of small-amplitude solutions? The an- 
swer is that gravity renders the zero-amplitude, i.e. pla- 
nar, solution more stable and hence larger amplitudes are 
needed for true structures to become stable. 

Hence, we conclude that in confined systems under 
gravity or a similar body force (it has been shown that 
in directional solidification a temperature gradient acts 
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just as a strong effective gravity field |21lf7] ) stable steady 
states may exist at large amplitude and be absent at small 
ones. 

Below the instability threshold, i.e., for parameters 
where the planar interface is stable, the system may be 
forced into a large amplitude state by a sufficiently strong 
perturbation. This clearly calls for numerical simulations 
and experimental attempts at creating these states. 

In extended systems, the absence of stable steady-state 
solutions at large wavelength as well as the numerical ev- 
idence from time-dependent simulations 10,14 suggest 
that the cusp singularity is indeed reached in finite time. 
This is of course a statement within linear elasticity the- 
ory. It means that stresses would increase beyond all limits 
in the minimimum of a groove, if linear elasticity held all 
the time. If linear elasticity were valid up to the fracture 
threshold, one might conclude from this result that the 
Grinfeld instability would inevitably lead to fracture in 
such a situation. However, the answer to this question is 
beyond the scope of this paper, as it is obvious that at 
sufficiently large stresses plasticity must be taken into ac- 
count. Phase-field simulations containing an inherent yield 
stress in the model jT3] suggest that indeed cracking is a 
likely scenario in sufficiently extended systems. 

This work was supported by the Deutsche Forsctiungsgemein- 
schaft under Grant No. Ka 672/4-2 and FOR 301/2-1, which 
is gratefully acknowledged. In addition, we acknowlede travel 
grants by PROCOPE, Grant No. 9619897 (DAAD, Germany) 
and 97176 (APAPE, France), enabling a closer collaboration 
between the two groups involved in this work. 

A Calculation of the strain energy density 
along a multi-cycloid surface 

Instead of calculating the strain energy density inside the 
bulk as in ^^li we will only describe here how to calculate 
this energy density w{s) at the surface. We carry out the 
calculation for general A^-cycloids. 

In general, the elastic energy density of a solid sub- 
mitted to plane strain can be written in terms of the two- 
dimensional stress tensor as 

1 + ^ ^ 2 \ /Ar\ 

W = [<^^J'^^j -''(^kk) > (45) 

where summation over repeated subscripts is implied. At 
the interface, aij is diagonal with elements au and (T,m- 
Since au = Tr cr — ann and because in our normalization 
Cnn = 0, the elastic energy density can be expressed by 
the trace of the stress tensor alone, which allows us to deal 
with a single scalar. Thus the strain energy density takes 
the form 

1 2 

w{x,y) ^ -{Tia) , (46) 
where Trcr is in our non-dimensional scalings equal to 1+ 

The idea is to employ a mapping of the half-plane 
bounded from above by our multi-cycloid onto the area 



below the real axis using the analytic function 

N 

z = a;(c) = c - ^ y ^e^^'^^ (47) 
^-^ nk 

n=l 

where <; = + irj. This defines a mapping of the domain 
'^{z) < C to the lower <j half plane, as can be seen easily 
by restricting to <j = ^ which shows that the interface 
is mapped to the real axis. In order to solve the elastic 
problem we have to satisfy ^1] 

MO + HO - 5111 +m) - ^^^^ , (48) 

^(0 2 

where 0o a-nd ipo are modified Goursat functions. These 
functions must be analytic functions of for 77 — > — cxd, so 
it has to be established that iPq{0 contains no exponen- 
tials increasing for rj —00 when ^ is replaced by ^ -|- irj. 
Since V'o(0 is the complex conjugate of a function that 
is analytic in the lower half plane, it must be analytic in 
the upper half plane, which means that terms of the form 
exp(m^) are allowed whereas exp(— m^) are not (for de- 
tails see 211). For brevity, we will designate the forbidden 
terms as "negative exponentials" . Technically, we make an 
ansatz for (/'0(C)- 

. N 

'^o(0 = ^E""e""'^' (49) 

n=l 

where we may assume q;„ to be real (which will be justified 
later). Now let us simplify Eq. (|48|l . We have 

HO - m] = 2z3 (c.(0) = E (50) 

and hence 

Via the choice of a„ we have to establish that the right 
hand side of Eq. (|51|l contains no negative exponentials. 
Let us further simplify the representation. We have 

N 

^(O^E"""^™'^'^ (52a) 

n=l 

N 

^(0 = 1 - E '^"e"'^ . (52b) 

n=l 

For the division, we use the common expression for the 
quotient of two series (^2;P- 28): Let 

51 = 1 4- aix + a2x'^ + a^x^ -I- . . . , 

52 ^ I + bix + h2X^ + 630:^ + ... , 

53 = 1 + CiX + C2X'^ + C^X^ + . . . , 

54 = 1 + dix + d2X^ + d^x^ + . . . , 

55 = 1 + eix + £22;^ + e.y,x^ + . . . , 
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with 
that is 



a„ = na. 
We then need to calculate 



X = e 



S3 



S2 



S5 = S3 - S4 = = . 

The coefficients c„ and (i„ are given by the recursion 
ci = ai - 61 , 



dn 



-hi , 



bn + bjdn^j 
J = l 



This leads to 



— ^ ^ bj 



(53) 



We then have, after introducing eo — 1/2, 



^ (U 



(54) 



n=0 



Therefore, 



Co 



"2 ' 



n-l 



n = 2...iV. 



(55) 



Negative exponentials on the right-hand side of (|51|l can 
only result from the negative exponentials in 



We cut out the relevant parts from Eq. (|F1T) , consisting of 
negative exponentials and require them to become zero, 
to compute the coefficients determining (f>Q in terms of the 

Pn 



N-l 




N 

E 



(56) 



with eo = s-nd the prefactor ik^^ dropped. Now we sort 
the terms in this equation by exponentials which finally 
gives us the system 



AT- 



Pn+j 



a„=^ej^^, n=l...N. 



n + j 



(57) 



Note that the contain the ai as well, so the equations 
are not a simple recursive scheme but a linear system of 
equations for the Q!„. Now we get back to Tr cr, which can 
be written as 

Tr..l + 4^iiK(f±4g^ (58) 

The denominator has already been written down during 
the calculation of the arc length. The numerator can be 
simplified in a similar manner. After all simplifications 
have been performed, we finally get 

/ N 1 . wi(s) 

Ms) = :j 1 + 4 ^ 

2 \ 1+ W2{S) 

N 

Wl{s) = Y (Cn - Pn) 

n=\ 

N-l N-n (59) 

~ E E (^j"j+«(-?' +r>')+ Pj+najj) 

ra=l j = l 

N N-l N-n 

W2{s) = Y Pn{Pn " 2C„) + 2 ^ C„ ^ P]P]+n ■ 
n— 1 n—1 j — 1 



B Details of the numerical solution of 

Eq. (I33D 

Let us first repeat the system of partial differential equa- 
tions igU). 



dPj 



(60) 



Herein, the energy changes entering the equation are inte- 
grals over certain rational terms containing trigonometric 
functions and the vector of amplitudes p = (pi . . . p^) as 
well as the wave number k and the gravity parameter I12. 

We first reformulate the problem H6U|I by considering 
the fact that k is contained in the energy terms as a pref- 
actor only (compare Eqs. (03) and gSJ)^. Then we 
get a simplified problem with modified E terms: 



d 
dpj 



\e,{p) + E,[p) + ^E^{p) 



= 0. 



(61) 



The question is now to specify which subset of the solu- 
tion manifold is required. We concentrate on fixed physical 



^ De facto we have also transformed the integration variable 
^ ^ X = fc^, changing the integration interval to [0, tt]. 
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system parameters, i.e., constant I12 in order to produce 
the lines shown in figure El 

Next a suitable numerical method has to be chosen. 
Solutions are known to satisfy p « for a ~ 0, and 
solutions starting a branch are given by linear stability 
analysis: fcstait = 1 + \/l - ^12; Pj, start = 0, j = 1 . . . TV. 

As solutions along a branch are expected to change 
continuously, we can implement a Newton Raphson al- 
gorithm and move along the selected branch by varying a 
parameter. Yet, already in the monocycloid model we find 
that some of the curves are multi-valued with respect to 
k. This renders it unfavourable to use fc as a fixed param- 
eter in that scheme and to solve for pi , because the exact 
turning points are unknowns. 

It turns out that for all branches situated between 
A; = 1 and k = 2, and with I12 & [0, 1], the curves behave 
monotonously as a function of pi up to the cusp. These 
are the cases exhibited in figure El Therefore, in order to 
solve the system of equations (|61|l for j = 1 ... TV, we keep 
pi fixed instead of k and take {k, p2 ■ . ■ pn} as our set of 
variables to be determined by the iteration. This results 
in a modified Jacobian containing terms {d^E /dpndk, 
d^E/dpndpm} {m = 2...N,n^l...N). 

As initial guess for the first non-zero solution of a 
branch, we choose the upper marginal wavenumber from 
linear theory for k, a small value of pi and set all other 
Pi equal to zero. After having found the first solution, we 
move along the solution branch towards larger pi values in 
steps of typically Api = 0.001. Consecutive data sets are 
estimated by forward differences using up to six solution 
points and then iterated until the 2-norm of the vector of 
the changes remains below a threshold of typically 10""'^°. 

It should be emphasized that the derivatives of Eqs. l|in|) . 
H42() and H43(l can be given analytically (they are omitted 
here because they are rather lengthy expressions), and so 
the Newton Raphson algorithm can be programmed with 
an (up to quadrature) exact Jacobi matrix. This makes 
the code converge extremely fast. 

A more thorough investigation of the solution mani- 
fold with respect to small wavenumbers and I12 > 1 goes 
beyond the scope of this paper because monotony consid- 
erations do not apply in this range and multi-valuedness 
of the manifold can appear which is connected to coars- 
ened solutions. These considerations will be presented in 
a separate work. 
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